library(fixest) # to demean variables
library(multiwayvcov) # to cluster SEs
library(sandwich) # to get HC2 SEs
library(stargazer) # to generate tables
library(mvtnorm) # to simulate betas and vcov matrix

rm(list = ls())

# load data
load(file = 'chDataWomenDistrict.RData')

# Put each dataset in a different object
ch_reform <- ch_list[[1]]
ch_refTerm <- ch_list[[2]]
ch_pre_post <- ch_list[[3]]
ch_refTermMChanges <- ch_list[[4]]

# "Gonzalo Fuenzalida" and 'Javier Macaya' changed district after the reform
ch_reform <- subset(ch_reform, !(bill_author %in% c("Gonzalo Fuenzalida", 'Javier Macaya')))
ch_refTerm <- subset(ch_refTerm, !(bill_author %in% c("Gonzalo Fuenzalida", 'Javier Macaya')))
ch_pre_post <- subset(ch_pre_post, !(bill_author %in% c("Gonzalo Fuenzalida", 'Javier Macaya')))
ch_refTermMChanges <- subset(ch_refTermMChanges, !(bill_author %in% c("Gonzalo Fuenzalida", 'Javier Macaya')))

# # Solo women
ch_reform <- subset(ch_reform, n_of_women == 1 | female == 0)
ch_refTerm <- subset(ch_refTerm, n_of_women == 1 | female == 0)
ch_pre_post <- subset(ch_pre_post, n_of_women == 1 | female == 0)
ch_refTermMChanges <- subset(ch_refTermMChanges, n_of_women == 1 | female == 0)

################################# 
####### Models using M ##########
#################################

# M1: G*M C1,C2,C3 (all data)
data_m1 <- fixest::demean(women_pct + magAnalysis + female + magAnalysisWoman
	~ session, data = ch_refTerm)
data_m1 <- cbind(data_m1, bill_author = ch_refTerm$bill_author)
m1 <- lm(women_pct ~ magAnalysis + female + magAnalysisWoman - 1, data = data_m1)
vcov1 <- cluster.vcov(m1, cluster = data_m1$bill_author)
s1  <- sqrt(diag(vcov1))

# M2: G*M C1,C2,C3 (M changes at the middle of C2)
data_m2 <- fixest::demean(women_pct + magAnalysis + female + magAnalysisWoman
	~ session + reform, data = ch_refTermMChanges)
data_m2 <- cbind(data_m2, bill_author = ch_refTermMChanges$bill_author)
m2 <- lm(women_pct ~ magAnalysis + female + magAnalysisWoman - 1, data = data_m2)
vcov2 <- cluster.vcov(m2, cluster = data_m2$bill_author)
s2  <- sqrt(diag(vcov2))

# M3: G*M C1,C3
data_m3 <- fixest::demean(women_pct + magAnalysis + female + magAnalysisWoman
	~ session, data = subset(ch_refTerm, session %in% c('2010-2014', '2018-2022')))
data_m3 <- cbind(data_m3, bill_author = subset(ch_refTerm, session %in% c('2010-2014', '2018-2022'))$bill_author)
m3 <- lm(women_pct ~ magAnalysis + female + magAnalysisWoman - 1, data = data_m3)
vcov3 <- cluster.vcov(m3, cluster = data_m3$bill_author)
s3  <- sqrt(diag(vcov3))

# M4: C2
inSample <- names(table(ch_pre_post$bill_author)[table(ch_pre_post$bill_author) == 2])
ch_pre_post <- subset(ch_pre_post, bill_author %in% inSample)
data_m4 <- fixest::demean(women_pct + magAnalysis + female + magAnalysisWoman
	~ bill_author, data = ch_pre_post)
data_m4 <- cbind(data_m4, bill_author = ch_pre_post$bill_author)
m4 <- lm(women_pct ~ magAnalysis + magAnalysisWoman - 1, data = data_m4)
vcov4 <- cluster.vcov(m4, cluster = data_m4$bill_author)
s4  <- sqrt(diag(vcov4))

# M6: G*M C1, C2, C3 (M changes at the middle of C2) (only repeated legislators)
data_m6 <- fixest::demean(women_pct + magAnalysis + magAnalysisWoman
	~ session + reform + bill_author, data = subset(ch_refTermMChanges, allSessions == 1))
data_m6 <- cbind(data_m6, bill_author = subset(ch_refTermMChanges, allSessions == 1)$bill_author)
m6 <- lm(women_pct ~ magAnalysis + magAnalysisWoman - 1, data = data_m6)
vcov6 <- cluster.vcov(m6, cluster = data_m6$bill_author)
s6  <- sqrt(diag(vcov6))


# M8: C3
data_m8 <- subset(ch_refTerm, session == '2018-2022')
m8 <- lm(women_pct ~ magAnalysis + female + magAnalysisWoman, data = data_m8)
vcov8 <- vcovHC(m8, 'HC2')
s8  <- sqrt(diag(vcov8))

# Generate Table I.33
stargazer(m1, m2, m3, m4, m6, m8, 
	se = list(s1, s2, s3, s4, s6, s8), 
	covariate.labels = c('M', 'Woman', 'M x Woman'), 
	add.lines = list(c('FE by Legislative Term', 'Yes', 'Yes', 'Yes', 'No', 'Yes', 'No'),
		c('FE by Reform', 'No', 'Yes', 'No', 'No', 'Yes', 'No'),
		c('FE by Legislator', 'No', 'No', 'No', 'Yes', 'Yes', 'No')),
	dep.var.labels = "Bills' Portfolio on Women's Issues (%)",
	omit.stat = c("ser",'f'),
	star.cutoffs = c(0.10, 0.05, 0.01),
	no.space = TRUE,
	single.row = FALSE,	
	notes.append = FALSE,
	title = "Association Between Legislative Portfolio, Gender, and District Magnitude--2014-2019  Cámara de Diputadas y Diputados de Chile (Only Woman)",
	label = 't:onlyWoman'
)

### Substantive Effect
## Scenarios: Women, M = 2
W2 <- data.frame(mag = 2, female = 1, magAnalysisWoman = 2 * 1)
## Scenarios: Women, M = 3 (post-reform)
W3 <- data.frame(mag = 3, female = 1, magAnalysisWoman = 3 * 1)
## Scenarios: Men, M = 2
M2 <- data.frame(mag = 2, female = 0, magAnalysisWoman = 2 * 0)
## Scenarios: Men, M = 3 (post-reform)
M3 <- data.frame(mag = 3, female = 0, magAnalysisWoman = 3 * 0)
# Women and Men Data
Wdata <- rbind(W2, W3)
Mdata <- rbind(M2, M3)

## Sims
set.seed(123)
sims1 <- rmvnorm(n = 1000, mean = coef(m1), sigma = vcov1)
sims2 <- rmvnorm(n = 1000, mean = coef(m2), sigma = vcov2)
sims3 <- rmvnorm(n = 1000, mean = coef(m3), sigma = vcov3)
sims4 <- rmvnorm(n = 1000, mean = coef(m4), sigma = vcov4)
sims6 <- rmvnorm(n = 1000, mean = coef(m6), sigma = vcov6)
sims8 <- rmvnorm(n = 1000, mean = coef(m8), sigma = vcov8)

## Calculate Predicted Values
# Women
pred1_W <- apply((sims1 %*% t(Wdata[2,])) - (sims1 %*% t(Wdata[1,])), 2, quantile, c(0.975, 0.5, 0.025))
pred2_W <- apply((sims2 %*% t(Wdata[2,])) - (sims2 %*% t(Wdata[1,])), 2, quantile, c(0.975, 0.5, 0.025))
pred3_W <- apply((sims3 %*% t(Wdata[2,])) - (sims3 %*% t(Wdata[1,])), 2, quantile, c(0.975, 0.5, 0.025))
pred4_W <- apply((sims4 %*% t(Wdata[2, -2])) - (sims4 %*% t(Wdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025))
pred6_W <- apply((sims6 %*% t(Wdata[2, -2])) - (sims6 %*% t(Wdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025))
pred8_W <- apply((sims8 %*% t(cbind(1, Wdata)[2,])) - (sims8 %*% t(cbind(1, Wdata)[1,])), 2, quantile, c(0.975, 0.5, 0.025))

# Men
pred1_M <- apply((sims1 %*% t(Mdata[2,])) - (sims1 %*% t(Mdata[1,])), 2, quantile, c(0.975, 0.5, 0.025))
pred2_M <- apply((sims2 %*% t(Mdata[2,])) - (sims2 %*% t(Mdata[1,])), 2, quantile, c(0.975, 0.5, 0.025))
pred3_M <- apply((sims3 %*% t(Mdata[2,])) - (sims3 %*% t(Mdata[1,])), 2, quantile, c(0.975, 0.5, 0.025))
pred4_M <- apply((sims4 %*% t(Mdata[2, -2])) - (sims4 %*% t(Mdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025))
pred6_M <- apply((sims6 %*% t(Mdata[2, -2])) - (sims6 %*% t(Mdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025))
pred8_M <- apply((sims8 %*% t(cbind(1, Mdata)[2,])) - (sims8 %*% t(cbind(1, Mdata)[1,])), 2, quantile, c(0.975, 0.5, 0.025))


# Building blocks of table 2
bill_hat_women <- c(round(pred1_W[2,], 2), round(pred2_W[2,], 2), round(pred3_W[2,], 2),
	round(pred4_W[2,], 2), round(pred6_W[2,], 2), round(pred8_W[2,], 2))
ci_women <- c(paste0('(', round(pred1_W[3,], 2), ', ', round(pred1_W[1,], 2), ')'), 
	paste0('(', round(pred2_W[3,], 2), ', ', round(pred2_W[1,], 2), ')'), 
	paste0('(', round(pred3_W[3,], 2), ', ', round(pred3_W[1,], 2), ')'), 
	paste0('(', round(pred4_W[3,], 2), ', ', round(pred4_W[1,], 2), ')'), 
	paste0('(', round(pred6_W[3,], 2), ', ', round(pred6_W[1,], 2), ')'), 
	paste0('(', round(pred8_W[3,], 2), ', ', round(pred8_W[1,], 2), ')'))


bill_hat_men <- c(round(pred1_M[2,], 2), round(pred2_M[2,], 2), round(pred3_M[2,], 2),
	round(pred4_M[2,], 2), round(pred6_M[2,], 2), round(pred8_M[2,], 2))
ci_men <- c(paste0('(', round(pred1_M[3,], 2), ', ', round(pred1_M[1,], 2), ')'), 
	paste0('(', round(pred2_M[3,], 2), ', ', round(pred2_M[1,], 2), ')'), 
	paste0('(', round(pred3_M[3,], 2), ', ', round(pred3_M[1,], 2), ')'), 
	paste0('(', round(pred4_M[3,], 2), ', ', round(pred4_M[1,], 2), ')'), 
	paste0('(', round(pred6_M[3,], 2), ', ', round(pred6_M[1,], 2), ')'), 
	paste0('(', round(pred8_M[3,], 2), ', ', round(pred8_M[1,], 2), ')'))

tableJ36 <- rbind("Women "= bill_hat_women, " " = ci_women, 
	"Men "= bill_hat_men, " " = ci_men)
colnames(tableJ36) <- paste0('Model ', c(1:4, 6, 8))
stargazer(tableJ36,
	title = "Predicted percentage point change in the cosponsorship portfolio on women’s issues when M increases by one")
